% Part 2: MCMC Methods
% This file is included by bayesian_vs_3dvar_cn_main.tex

\section{马尔可夫链蒙特卡罗方法}

\subsection{MCMC 基本原理}
当直接求解析解计算量过大或需要完整后验分布时,可采用 MCMC 方法从后验分布中采样。MCMC 通过构造一个平稳分布为目标后验分布的马尔可夫链,经过足够长的迭代后,链的状态收敛到目标分布。

{\color{red}\textbf{[采样 vs 优化]：MCMC 与梯度优化代表两种不同的哲学。优化寻找后验众数(单点估计)，采样探索整个后验分布(概率估计)。MCMC 的优势在于：(1)无需计算梯度；(2)可处理非高斯、多峰分布；(3)自然量化不确定性。代价是需要大量样本($N \sim 10^4$-$10^6$)才能获得可靠估计。}}

\subsection{Metropolis-Hastings 算法}
Metropolis-Hastings (MH) 是最基本的 MCMC 算法。在第 $t$ 步,给定当前状态 $\mathbf{x}^{(t)}$:
\begin{enumerate}
\item 从提议分布 $q(\mathbf{x}' | \mathbf{x}^{(t)})$ 中生成候选状态 $\mathbf{x}'$;
\item 计算接受概率
\begin{equation}
\alpha = \min\left(1, \frac{p(\mathbf{x}' | \mathbf{y}^{\text{obs}}) q(\mathbf{x}^{(t)} | \mathbf{x}')}{p(\mathbf{x}^{(t)} | \mathbf{y}^{\text{obs}}) q(\mathbf{x}' | \mathbf{x}^{(t)})}\right);
\end{equation}
\item 以概率 $\alpha$ 接受 $\mathbf{x}^{(t+1)} = \mathbf{x}'$,否则保持 $\mathbf{x}^{(t+1)} = \mathbf{x}^{(t)}$。
\end{enumerate}

对于对称提议分布 $q(\mathbf{x}' | \mathbf{x}) = q(\mathbf{x} | \mathbf{x}')$,接受概率简化为
\begin{equation}
\alpha = \min\left(1, \frac{p(\mathbf{x}' | \mathbf{y}^{\text{obs}})}{p(\mathbf{x}^{(t)} | \mathbf{y}^{\text{obs}})}\right).
\end{equation}

{\color{red}\textbf{[接受概率的智慧]：MH 接受准则的精妙之处在于：(1)总是接受``更好''的状态($p(\mathbf{x}') > p(\mathbf{x}^{(t)})$)；(2)以概率 $p(\mathbf{x}')/p(\mathbf{x}^{(t)})$ 接受``更差''的状态。这允许链逃离局部极值，探索整个状态空间。关键是只需计算概率比值，无需归一化常数！}}

\subsection{随机游走 Metropolis}
最简单的提议为 $\chi$ 空间随机游走:
\begin{equation}
\chi' = \chi^{(t)} + \sigma \boldsymbol{\epsilon}, \quad \boldsymbol{\epsilon} \sim N(0, \mathbf{I}),
\end{equation}
其中 $\sigma$ 为步长参数。最优接受率约为 23\%,步长需满足 $\sigma^2 \propto d^{-1}$($d$ 为维度)。

{\color{red}\textbf{[维度诅咒]：随机游走的效率随维度急剧下降。步长必须按 $d^{-1}$ 缩放以维持合理接受率，导致混合时间按 $d^2$ 增长。对于 $d \sim 10^4$ 的大气反演问题，随机游走几乎不可用，需要更高级的方法如 pCN 或 HMC。}}

\subsection{预条件 Crank-Nicolson 算法}
利用先验协方差结构的 pCN 算法在 $\chi$ 空间提议
\begin{equation}
\chi' = \sqrt{1-\beta^2}\chi^{(t)} + \beta\boldsymbol{\epsilon}, \quad \boldsymbol{\epsilon} \sim N(0, \mathbf{I}),
\end{equation}
其中 $\beta \in (0,1)$ 控制步长。pCN 的关键优势是保持先验分布不变,接受率与维度无关,适用于高维问题($n > 10^4$)。

{\color{red}\textbf{[维度独立性的奇迹]：pCN 通过在先验支撑集内提议，避免了维度诅咒。接受率仅依赖于似然比，与维度 $d$ 无关！这源于提议保持先验不变的性质：$\chi' \sim N(0,\mathbf{I})$ 与 $\chi \sim N(0,\mathbf{I})$ 同分布。这使 pCN 成为高维($d > 10^4$)大气反演的首选算法。}}

\subsection{Metropolis 调整 Langevin 算法}
MALA 利用梯度信息引导采样方向:
\begin{equation}
\chi' = \chi^{(t)} + \frac{\sigma^2}{2}\nabla_\chi \log p(\chi^{(t)} | \mathbf{y}^{\text{obs}}) + \sigma\boldsymbol{\epsilon},
\end{equation}
其中 $\chi$ 空间的梯度为
\begin{equation}
\nabla_\chi \log p(\chi | \mathbf{y}^{\text{obs}}) = \mathbf{B}^{1/2}\mathbf{H}^T\mathbf{R}^{-1}(\mathbf{y}^{\text{obs}} - \mathbf{H}\mathbf{x}) - \chi.
\end{equation}
MALA 的收敛速度优于随机游走,步长缩放关系为 $\sigma^2 \propto d^{-1/3}$。

{\color{red}\textbf{[梯度引导的力量]：MALA 通过梯度``指向''高概率区域，比盲目随机游走高效得多。其提议是 Langevin 扩散的离散化，自然平衡漂移(梯度)与扩散(随机噪声)。相比随机游走的 $d^{-1}$ 缩放，MALA 的 $d^{-1/3}$ 缩放显著减轻维度依赖，但仍不如 pCN 的维度独立性。}}

\subsection{收敛诊断}
\subsubsection{Gelman-Rubin 统计量}
对于 $m$ 条独立链,每条长度 $n$,计算链间方差 $B$ 和链内方差 $W$:
\begin{equation}
B = \frac{n}{m-1}\sum_{j=1}^m(\bar{\theta}_j - \bar{\bar{\theta}})^2, \quad W = \frac{1}{m}\sum_{j=1}^m s_j^2,
\end{equation}
其中 $\bar{\theta}_j$ 为第 $j$ 条链的均值,$s_j^2$ 为方差。$\hat{R}$ 统计量定义为
\begin{equation}
\hat{R} = \sqrt{\frac{n-1}{n} + \frac{1}{n}\frac{B}{W}}.
\end{equation}
当 $\hat{R} < 1.1$ 时认为已收敛,理想情况下 $\hat{R} \approx 1$。

{\color{red}\textbf{[收敛判据]：$\hat{R}$ 比较``链间变异''与``链内变异''。若链未收敛，不同初值的链会探索不同区域，导致 $B \gg W$ 和 $\hat{R} \gg 1$。收敛后，所有链混合在同一分布中，$B \approx W$ 和 $\hat{R} \approx 1$。阈值 $\hat{R} < 1.1$ 是经验标准，严格应用建议 $\hat{R} < 1.01$。}}

\subsubsection{有效样本量}
考虑自相关后的有效样本量为
\begin{equation}
\text{ESS} = \frac{N}{1 + 2\sum_{k=1}^K \rho_k},
\end{equation}
其中 $\rho_k$ 为滞后 $k$ 的自相关系数,$K$ 为截断点。通常要求 ESS $> 400$ 以获得可靠估计。

{\color{red}\textbf{[自相关的代价]：MCMC 样本是相关的，不同于独立同分布(i.i.d.)样本。强自相关($\rho_k$ 衰减慢)意味着连续样本几乎相同，信息冗余。ESS 量化``等效独立样本数''：若 $N=10^6$ 但 ESS$=400$，则仅相当于400个独立样本的统计效力。提高 ESS 需减小自相关(更好的提议)或增加总样本数。}}

\subsection{后验分布估计}
经过 $N$ 次迭代后,丢弃前 $N_{\text{burn}}$ 个老化样本,后验均值和协方差估计为
\begin{equation}
\mathbb{E}[\mathbf{x} | \mathbf{y}^{\text{obs}}] \approx \frac{1}{N - N_{\text{burn}}}\sum_{t=N_{\text{burn}}+1}^N \mathbf{x}^{(t)},
\end{equation}
\begin{equation}
\text{Cov}[\mathbf{x} | \mathbf{y}^{\text{obs}}] \approx \frac{1}{N - N_{\text{burn}} - 1}\sum_{t=N_{\text{burn}}+1}^N (\mathbf{x}^{(t)} - \bar{\mathbf{x}})(\mathbf{x}^{(t)} - \bar{\mathbf{x}})^T.
\end{equation}

